################################################################################
# Stat-Tutorial-02-Distributions.r                                             #
#                                                                              #
# This is a simple tutorial to introduce basic calculations related to key     #
# probability distributions using R.                                           #
#                                                                              #
# R. Labouriau                                                                 #
# Last revised: Fall 2020                                                      #
################################################################################
# Copyright © 2019 by Rodrigo Labouriau

################################################################################
# Binomial distribution
################################################################################

# Returns the probability function of a binomial distribution
# with 10 independent trials (i.e. size = 10) each with probability 0.5 (prob = 0.5)
# of success. Answers to the question: What is the probability
# of tossing a fair coin 10 times and get only 1 head.
dbinom(x=1, size=10, prob=0.5) 

# Returns the c.d.f. of a binomial distribution
# with 10 independent trials (i.e. size = 10) each with probability 0.5 (prob = 0.5)
# of success. That is, if X ~ Bi (10, 0.5), the function returns 
# the P(X lexx or equal 8). Answers to the question: What is the probability
# of tossing a fair coin 10 times and get less or equal 8 heads.
pbinom(q=8, size=10, prob=0.5)

# Returns x such that P(X<=x) = 0.1, if X ~ Bi (10, 0.5)
qbinom(p=0.1, size=10, prob=0.5)  

# Generates 5 "random" values drawn fron a Bi (10, 0.5) distribution
rbinom(n=5, size=10, prob=0.5)

# Here, we make some illustrative figures
B <- rbinom(n=1000, size=10, prob=0.5)
hist(B, col="lightblue", freq=F, ylim=c(0,1/2), main="")
points(0:10, dbinom(x=0:10, size=10, prob=0.5), type="h", lwd=5)

B <- rbinom(n=100, size=10, prob=0.5)
hist(B, col="lightblue", freq=F,)
points(0:10, dbinom(x=0:10, size=10, prob=0.5), type="h", lwd=5)

################################################################################
# QUESTIONS:
# 1) What is the probability of tossing a fair coin 100 times and getting 
# exactly 30 heads? 
#
# 2) What is the probability of tossing a fair coin 100 times and getting 
# less than 30 heads?
#
# 3) What is the probability of tossing a fair coin 100 times and getting 
# less than 30 tails?
#
# Optional additional question:
# 4*) Suppose that you toss a coin 10 times and get 9 tails. What is the 
# probability that this happens if the coin is fair? What is the probability
# that this happens if the coin is such that the probability of tail is 3/4?
# Or if the probability of tail is 0.9? Or 0.91? Or 0.92 ?
#
################################################################################

################################################################################
# Normal distributions                                                         #
################################################################################

# Simulating a (pseudo) random variable drawn from a normal distribution with
# mean 0 and standard error 1 and 10 values
rnorm(n=10, mean=0, sd=1)

rnorm(n=10, mean=0, sd=1)

rnorm(10, 0, 1)

Z <- rnorm(n=100, mean=0, sd=1)
plot(Z)

plot(Z, type="l")

Z <- rnorm(n=100, mean=0, sd=1)
plot(Z, type="l")

hist(Z, col="lightblue")

pnorm(1.5, 1, 2) # Returns the value of the normal(1,2) cdf at x=1.5

qnorm(0.2, 1, 2) # Returns x such that P(X<=x) = 0.2 where X is a
# normal(1,2) random variable (i.e. returns the cdf)

dnorm (0, 1, 2)  #  Returns the value of the density of a normal distribution
#  with mean 1 and sd 2 evaluated at the point 0

# We draw an illustrative figure
Z <- rnorm(n=10000, mean=0, sd=1)

hist(Z, col="lightblue", freq=F, main="10000 simulated standard normal observations")
curve(dnorm(x, mean=0,sd=1), from=-4, to=4, add=T, lwd=3)

################################################################################
# Uniform distribution
################################################################################

dunif(x=0.3, min=0, max=1) # Returns the density of a uniform distribution at 0.3
punif(0.3, min=0, max=1)   # Returns the c.d.f. of a uniform distribution at 0.3
qunif(p=0.1, min=0, max=1) # returns the quantiles (p=0.1) of an uniform distr.
runif(n=5, min=0, max=1)   # Simulates 5 values drawn from a uniform distribution

hist(runif(1000), col="lightblue", main="")

hist(runif(10), col="lightblue", main="")

################################################################################
# Poisson distribution
################################################################################

# Returns the density of a Poisson distribution with lambda=5 at 1
dpois(1, lambda=5)  
# Returns the c.d.f. of a Poisson distribution with lambda=5 at 1
ppois(1, lambda=5)
# Returns the p-quantile of a Poisson distribution with lambda=5 
qpois(p, lambda)
# Generates 100 "random" values drawn from a Poisson distribution with
# lambda=10
rpois(n=100, lambda=10)

hist(rpois(n=1000, lambda=10), col="lightblue", main="", freq=F)
points(0:20, dpois(0:20, lambda=10), type="h", lwd=5)

Y <- rpois(n=1000, lambda=5)
hist(Y, col="lightblue", main="Wrong distribution", freq=F)
points(0:20, dpois(0:20, lambda=10), type="h", lwd=5)

hist(Y, col="lightblue", main="Right distribution", freq=F)
points(0:20, dpois(0:20, lambda=5), type="h", lwd=5)

# THE END!